home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / delaunay / quadedge.C < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  6.8 KB  |  246 lines
  1. #include <quadedge.h>
  2.  
  3. /*********************** Basic Topological Operators ************************/
  4.  
  5. Edge* MakeEdge()
  6. {
  7.     QuadEdge *ql = new QuadEdge;
  8.     return ql->e;
  9. }
  10.  
  11. void Splice(Edge* a, Edge* b)
  12. // This operator affects the two edge rings around the origins of a and b,
  13. // and, independently, the two edge rings around the left faces of a and b.
  14. // In each case, (i) if the two rings are distinct, Splice will combine
  15. // them into one; (ii) if the two are the same ring, Splice will break it
  16. // into two separate pieces.
  17. // Thus, Splice can be used both to attach the two edges together, and
  18. // to break them apart. See Guibas and Stolfi (1985) p.96 for more details
  19. // and illustrations.
  20. {
  21.     Edge* alpha = a->Onext()->Rot();
  22.     Edge* beta  = b->Onext()->Rot();
  23.  
  24.     Edge* t1 = b->Onext();
  25.     Edge* t2 = a->Onext();
  26.     Edge* t3 = beta->Onext();
  27.     Edge* t4 = alpha->Onext();
  28.  
  29.     a->next = t1;
  30.     b->next = t2;
  31.     alpha->next = t3;
  32.     beta->next = t4;
  33. }
  34.  
  35. void DeleteEdge(Edge* e)
  36. {
  37.     Splice(e, e->Oprev());
  38.     Splice(e->Sym(), e->Sym()->Oprev());
  39.     delete e->Qedge();
  40. }
  41.  
  42. /************* Topological Operations for Delaunay Diagrams *****************/
  43.  
  44. Subdivision::Subdivision(const Point2d& a, const Point2d& b, const Point2d& c)
  45. // Initialize a subdivision to the triangle defined by the points a, b, c.
  46. {
  47.     Point2d *da, *db, *dc;
  48.     da = new Point2d(a), db = new Point2d(b), dc = new Point2d(c);
  49.     Edge* ea = MakeEdge();
  50.     ea->EndPoints(da, db);
  51.     Edge* eb = MakeEdge();
  52.     Splice(ea->Sym(), eb);
  53.     eb->EndPoints(db, dc);
  54.     Edge* ec = MakeEdge();
  55.     Splice(eb->Sym(), ec);
  56.     ec->EndPoints(dc, da);
  57.     Splice(ec->Sym(), ea);
  58.     startingEdge = ea;
  59. }
  60.  
  61. Edge* Connect(Edge* a, Edge* b)
  62. // Add a new edge e connecting the destination of a to the
  63. // origin of b, in such a way that all three have the same
  64. // left face after the connection is complete.
  65. // Additionally, the data pointers of the new edge are set.
  66. {
  67.     Edge* e = MakeEdge();
  68.     Splice(e, a->Lnext());
  69.     Splice(e->Sym(), b);
  70.     e->EndPoints(a->Dest(), b->Org());
  71.     return e;
  72. }
  73.  
  74. void Swap(Edge* e)
  75. // Essentially turns edge e counterclockwise inside its enclosing
  76. // quadrilateral. The data pointers are modified accordingly.
  77. {
  78.     Edge* a = e->Oprev();
  79.     Edge* b = e->Sym()->Oprev();
  80.     Splice(e, a);
  81.     Splice(e->Sym(), b);
  82.     Splice(e, a->Lnext());
  83.     Splice(e->Sym(), b->Lnext());
  84.     e->EndPoints(a->Dest(), b->Dest());
  85. }
  86.  
  87. /*************** Geometric Predicates for Delaunay Diagrams *****************/
  88.  
  89. inline Real TriArea(const Point2d& a, const Point2d& b, const Point2d& c)
  90. // Returns twice the area of the oriented triangle (a, b, c), i.e., the
  91. // area is positive if the triangle is oriented counterclockwise.
  92. {
  93.     return (b.x - a.x)*(c.y - a.y) - (b.y - a.y)*(c.x - a.x);
  94. }
  95.  
  96. int InCircle(const Point2d& a, const Point2d& b,
  97.              const Point2d& c, const Point2d& d)
  98. // Returns TRUE if the point d is inside the circle defined by the
  99. // points a, b, c. See Guibas and Stolfi (1985) p.107.
  100. {
  101.     return (a.x*a.x + a.y*a.y) * TriArea(b, c, d) -
  102.            (b.x*b.x + b.y*b.y) * TriArea(a, c, d) +
  103.            (c.x*c.x + c.y*c.y) * TriArea(a, b, d) -
  104.            (d.x*d.x + d.y*d.y) * TriArea(a, b, c) > 0;
  105. }
  106.  
  107. int ccw(const Point2d& a, const Point2d& b, const Point2d& c)
  108. // Returns TRUE if the points a, b, c are in a counterclockwise order
  109. {
  110.     return (TriArea(a, b, c) > 0);
  111. }
  112.  
  113. int RightOf(const Point2d& x, Edge* e)
  114. {
  115.     return ccw(x, e->Dest2d(), e->Org2d());
  116. }
  117.  
  118. int LeftOf(const Point2d& x, Edge* e)
  119. {
  120.     return ccw(x, e->Org2d(), e->Dest2d());
  121. }
  122.  
  123. int OnEdge(const Point2d& x, Edge* e)
  124. // A predicate that determines if the point x is on the edge e.
  125. // The point is considered on if it is in the EPS-neighborhood
  126. // of the edge.
  127. {
  128.     Real t1, t2, t3;
  129.     t1 = (x - e->Org2d()).norm();
  130.     t2 = (x - e->Dest2d()).norm();
  131.     if (t1 < EPS || t2 < EPS)
  132.         return TRUE;
  133.     t3 = (e->Org2d() - e->Dest2d()).norm();
  134.     if (t1 > t3 || t2 > t3)
  135.         return FALSE;
  136.     Line line(e->Org2d(), e->Dest2d());
  137.     return (fabs(line.eval(x)) < EPS);
  138. }
  139.  
  140. /************* An Incremental Algorithm for the Construction of *************/
  141. /************************ Delaunay Diagrams *********************************/
  142.  
  143. Edge* Subdivision::Locate(const Point2d& x)
  144. // Returns an edge e, s.t. either x is on e, or e is an edge of
  145. // a triangle containing x. The search starts from startingEdge
  146. // and proceeds in the general direction of x. Based on the
  147. // pseudocode in Guibas and Stolfi (1985) p.121.
  148. {
  149.     Edge* e = startingEdge;
  150.  
  151.     while (TRUE) {
  152.         if (x == e->Org2d() || x == e->Dest2d())
  153.             return e;
  154.         else if (RightOf(x, e))
  155.              e = e->Sym();
  156.         else if (!RightOf(x, e->Onext()))
  157.              e = e->Onext();
  158.         else if (!RightOf(x, e->Dprev()))
  159.              e = e->Dprev();
  160.         else
  161.             return e;
  162.     }
  163. }
  164.  
  165. void Subdivision::InsertSite(const Point2d& x)
  166. // Inserts a new point into a subdivision representing a Delaunay
  167. // triangulation, and fixes the affected edges so that the result
  168. // is still a Delaunay triangulation. This is based on the
  169. // pseudocode from Guibas and Stolfi (1985) p.120, with slight
  170. // modifications and a bug fix.
  171. {
  172.     Edge* e = Locate(x);
  173.     if ((x == e->Org2d()) || (x == e->Dest2d()))  // point is already in
  174.         return;
  175.     else if (OnEdge(x, e)) {
  176.         e = e->Oprev();
  177.         DeleteEdge(e->Onext());
  178.     }
  179.  
  180.     // Connect the new point to the vertices of the containing
  181.     // triangle (or quadrilateral, if the new point fell on an
  182.     // existing edge.)
  183.     Edge* base = MakeEdge();
  184.     base->EndPoints(e->Org(), new Point2d(x));
  185.     Splice(base, e);
  186.     startingEdge = base;
  187.     do {
  188.         base = Connect(e, base->Sym());
  189.         e = base->Oprev();
  190.     } while (e->Lnext() != startingEdge);
  191.  
  192.     // Examine suspect edges to ensure that the Delaunay condition
  193.     // is satisfied.
  194.     do {
  195.         Edge* t = e->Oprev();
  196.         if (RightOf(t->Dest2d(), e) &&
  197.             InCircle(e->Org2d(), t->Dest2d(), e->Dest2d(), x)) {
  198.                 Swap(e);
  199.                 e = e->Oprev();
  200.         }
  201.         else if (e->Onext() == startingEdge)  // no more suspect edges
  202.             return;
  203.         else  // pop a suspect edge
  204.             e = e->Onext()->Lprev();
  205.     } while (TRUE);
  206. }
  207.  
  208. /*****************************************************************************/
  209.  
  210. #include <gl.h>
  211. static unsigned int timestamp = 0;
  212.  
  213. void Subdivision::Draw()
  214. {
  215.     if (++timestamp == 0)
  216.         timestamp = 1;
  217.     startingEdge->Draw(timestamp);
  218. }
  219.  
  220. void Edge::Draw(unsigned int stamp)
  221. // This is a recursive drawing routine that uses time stamps to
  222. // determine if the edge has already been drawn. This is given
  223. // here for testing purposes only: it is not efficient, and for
  224. // large triangulations the stack might overflow. A better way
  225. // of doing this (and other traversals of the edges) is to maintain
  226. // a list of edges in the corresponding Subdivision object. This
  227. // list should be updated every time an edge is created or destroyed.
  228. {
  229.     if (Qedge()->TimeStamp(stamp)) {
  230.  
  231.         // Draw the edge
  232.         Point2d a = Org2d();
  233.         Point2d b = Dest2d();
  234.         bgnline();
  235.         v2d((double*)&a);
  236.         v2d((double*)&b);
  237.         endline();
  238.  
  239.         // visit neighbors
  240.         Onext()->Draw(stamp);
  241.         Oprev()->Draw(stamp);
  242.         Dnext()->Draw(stamp);
  243.         Dprev()->Draw(stamp);
  244.     }
  245. }
  246.